Preamble

This notebook deals with compiling all available estimates and measurements of vertical nitrate fluxes (in the Arctic Ocean).

All figures exported from this notebook are prefixed by FIGURE_NO3-COMP_.

In [1]:
%load_ext autoreload
%autoreload 2
%run imports.py

dims = dict(
    FN = hv.Dimension('FN', label='Nitrate flux', unit='mmol N m⁻² d⁻¹', range=(.005, 10)),
    no3_sfc_winter_uM = hv.Dimension('no3_sfc_winter_uM', label='Pre-bloom sfc. NO₃ ', unit='µM', range=(-.5, 13)),
    strat=hv.Dimension('strat', label='Brunt-Väisälä frequency', unit='s⁻²')
)

Maps

In [16]:
# df = gpd.read_file('../data/fn-compilation-database/fn-compilation.shp')
# df = gpd.read_file('../data/fn-compilation')
df = pd.read_csv('../data/fn-compilation.csv')
# df['geometry'] = df.centroid

# df.Longitude = df.geometry.x
# df.Latitude = df.geometry.y

# df = pd.DataFrame(df.drop(columns='geometry'))
# df['FN'] = df.FN.replace(0, 1e-10)

options_bk = [
    opts.Points(
        backend='bokeh',
        color=dim('FN'),
        clipping_colors={'min': 'grey', 'NaN': 'grey'},
        #hooks=[logcolor_ticks([1e-2, 1e-1, 1, 5])],
        tools=['hover'], cmap=default_cmap, colorbar=True, 
        size=dim('samplesize').categorize({'single': 9, 'aggregate': 20}),
        line_color='k',
        width=500, height=400, show_legend=False),
    opts.Labels(xoffset=5e3, yoffset=5e3),
    opts.Feature(scale='50m'),
]
             

options_mpl = translate_options(options_bk, bokeh2mpl)
options = options_mpl + options_bk

data = gv.Points(df, 
                 kdims=['Longitude', 'Latitude'], 
                 vdims=['FN','samplesize', 'Reference', 'Season'], 
                 crs=ccrs.PlateCarree()
                )

# coastline = gf.coastline.clone(extents=(-180, 70, 180, 90))

l = land * data #* gv.Labels(data)
# l = data
l = l.redim.range(Latitude=(60,90), Longitude=(-180, 180))

ll = (
    l
    .opts(*(opts_map+options), clone=True)
    .opts(
        opts.Points(
            backend='bokeh', 
            color=np.log10(dim('FN')).clip(min=-2), hooks=[logcolor_ticks([1e-2, 1e-1, 1, 5])]
        )
    )
)
hv.renderer('bokeh').theme = theme
hv.save(ll, '../nb_fig/FIGURE_FN-COMP_map.html')
ll
/Users/doppler/anaconda3/envs/nitrateflux/lib/python3.7/site-packages/holoviews/util/transform.py:334: RuntimeWarning: divide by zero encountered in log10
  data = o['fn'](*args, **kwargs)
/Users/doppler/anaconda3/envs/nitrateflux/lib/python3.7/site-packages/holoviews/util/transform.py:334: RuntimeWarning: divide by zero encountered in log10
  data = o['fn'](*args, **kwargs)
/Users/doppler/anaconda3/envs/nitrateflux/lib/python3.7/site-packages/holoviews/util/transform.py:334: RuntimeWarning: divide by zero encountered in log10
  data = o['fn'](*args, **kwargs)
/Users/doppler/anaconda3/envs/nitrateflux/lib/python3.7/site-packages/holoviews/util/transform.py:334: RuntimeWarning: divide by zero encountered in log10
  data = o['fn'](*args, **kwargs)
Out[16]:
In [3]:
def finalize(fig):
    ax = fig.axes[0]
    ax.background_img()
    
fig = mplrender_map(
    l
    .opts(
        *opts_map_mpl, *options_mpl,
        opts.Points(
            logz=True, 
            cmap=cmocean.cm.solar, zlim=(1e-2, 5), clipping_colors={'min': 'grey', 'NaN': 'grey'}, 
            backend='matplotlib'
        ),
        opts.Feature(scale='10m'),
        clone=True
    ),
    '../nb_fig/FIGURE_FN-COMP_map',
    hooks=[finalize],
)
fig
Out[3]:

GridMatrix (faceted)

def subset(N, season): df2=df.loc[(df.samplesize==N)&(df.Season==season)] if len(df2)>0: return gv.Points(df2, ['Longitude', 'Latitude'], ['FN', 'logFN','samplesize'], crs=ccrs.PlateCarree()) else: return gv.Points(None) griddict = {(season, N): subset(N, season) for N in ['single', 'aggregate'] for season in df.Season.unique()} fn = hv.HoloMap(griddict, ['season', 'samplesize']) fn_all = gv.Points(df, ['Longitude', 'Latitude'], crs=ccrs.PlateCarree(), group='all') l = geography * fn_all * fn l.opts(*options) l = l.opts(opts.Points('all', color='k', alpha=.4))l = geography * fn[['Perennial', 'Winter'], :].overlay().opts(title_format='Nitrate fluxes in winter') hv.save(l.opts(toolbar=None, clone=True), '../nb_fig/FIGURE_FN-COMP_map_winter.png') hv.save(l, '../nb_fig/FIGURE_FN-COMP_map_winter.html') lp = hv.render(l)

Seasonal cycle

In [4]:
df = pd.read_csv('../data/fn-compilation.csv')

df_nonwrap = df.loc[df.validity_start<df.validity_end]
df_wrap = df.loc[df.validity_start>df.validity_end]

df = pd.concat([
    df_nonwrap,
    df_wrap.assign(measured_end=13, validity_end=13),
    df_wrap.assign(measured_start=1, validity_start=1)
])

# reshape the start and end times of the different periods into a tidy tabular format 
# for easier ingestion in holoviews
def melt_df(x):
    dfr = df.melt(id_vars=['Reference', 'FN', 'samplesize', 'overturning'], value_vars=['measured_'+x, 'validity_'+x], 
                  value_name=x, var_name='type')
    dfr.loc[dfr.type.str.contains('measured'), 'type']='measured'
    dfr.loc[dfr.type.str.contains('valid'), 'type']='valid'
    return dfr

df = melt_df('start').join(melt_df('end')[['end']])

df = df.replace(dict(
    aggregate='Larger-scale aggregate', single='Few stations', overturning='Winter overturning', perennial='Perennial stratification', 
    measured='Measured', valid='Stipulated'
))

df.loc[df.start==df.end, 'end'] += 1

df.head()
Out[4]:
Reference FN samplesize overturning type start end
0 randelhoff2016vertical 0.300 Larger-scale aggregate Winter overturning Measured 5.0 9.0
1 randelhoff2016vertical 0.700 Larger-scale aggregate Winter overturning Measured 5.0 9.0
2 bourgault2011turbulent 0.120 Larger-scale aggregate Perennial stratification Measured 11.0 12.5
3 randelhoff2016regional 0.010 Larger-scale aggregate Perennial stratification Measured 4.0 5.0
4 randelhoff2016regional 0.015 Larger-scale aggregate Perennial stratification Measured 4.0 5.0
In [5]:
t=pd.date_range('2000-1-1', '2001-1-1', freq='2MS')
xticks = list(zip(t.month+.5, t.month_name().str[:3]))

od = {
    'overturning': {'color': ['orange', 'blue']},
    'type': {'line_width': [6, 3]},
    'samplesize': {'alpha': [.3, 1]}
}

options = (
    opts.NdOverlay(legend_position='right', show_legend=False, width=800),
    opts.Segments(logy=True, padding=.1, xticks=xticks, show_title=False, title_format='',
                  toolbar=None, yformatter='%g', **faceted_legend_opts(od), xlabel=''),
)

l = hv.Dataset(df.assign(fn2=df.FN), ['start', 'FN', 'end', 'fn2'] + ['overturning', 'type', 'samplesize']).to(hvu.Segments).overlay()
l = l.redim(**dims)
l = l.opts(*options) + faceted_legend(l, od, hw={'height':200, 'width':250})

fname = '../nb_fig/FIGURE_FN-COMP_chart_seasonal_cycle'
ll = l.opts(toolbar=None, clone=True)
hv.save(ll, fname, fmt='png')
save_bokeh_svg_multipanel(ll, fname+'.svg')
hv.save(l, fname, fmt='html')
l
Out[5]:

Surface nitratre concentration vs. stratification vs. bathymetry

def adjust_legend(fig): ax = fig.axes[0] leg = ax.get_legend() leg.set_title(None)

FN- winter NO3

In [6]:
df = pd.read_csv('../data/fn-compilation.csv')
df.overturning = df.overturning.str.capitalize()

options = [
    opts.Scatter(
        color='black',
        marker=hv.Cycle(['triangle', 'square']),
        show_grid=False, 
        legend_position='top_left', show_legend=True,
        height=350, width=350, #color=dim('label'),
        size=10, 
        logx=True, xformatter='%g'),
]

l = (
    hv.Dataset(df, ['FN', 'overturning'], 'no3_sfc_winter_uM')
    .to(hv.Scatter).overlay('overturning')
    .redim(**dims)
)

l = l.opts(*options).redim.range(FN=((5e-3,6))).relabel('A')

hv.renderer('bokeh').theme = theme
fname = '../nb_fig/FIGURE_FN-COMP_chart_vs_NO3_panelA'
hv.save(l.opts(toolbar=None, clone=True), fname, fmt='html')
hv.save(l.opts(toolbar=None, clone=True), fname, fmt='png')
save_bokeh_svg(l, fname+'.svg')
panel_A = l.clone()
panel_A
Out[6]:

Panel B: Surface nitratre vs. stratification grouped by bathymetry

In [7]:
unicode_dict = {i:j for i, j in zip(range(10), list('⁰¹²³⁴⁵⁶⁷⁸⁹'))}

def format_exponent_unicode(n):
    if n>0:
        return '10'+unicode_dict[n]
    elif n<0:
        return '10⁻'+unicode_dict[abs(n)]
In [8]:
df = (
    pd.read_csv('../data/no3-compilation/database-per-stn.csv', parse_dates=['date'])
    .rename(columns=dict(N2_30_100='strat'))
)
df = df.assign(logstrat=np.log10(df.strat))

df = df.loc[df.month.between(3, 5)]

ds = xr.load_dataarray('/Users/doppler/database/IBCAO/IBCAO_1min_bathy.nc')
df['btm'] = - ds.sel(lat=df.lat.values, lon=df.lon.values, method='nearest').values.diagonal()

# df = df.groupby('reg_name', as_index=False).mean()
df = df.dropna(subset=['reg_name'])
df['btm_cat'] = pd.cut(df.btm, bins=[0, 200, 1500, np.inf], labels=['Shelf', 'Slope', 'Basin'])

strat_bins = np.arange(-6, -3, .5)
df['strat_cat'] = pd.cut(df.logstrat, bins=strat_bins)

btm_cats = sc.data.keys()
# ds = hv.Dataset('logstrat', ['ntr0', 'btm_cat'])

def binavg(el):
    el = bin_average(el, bins=strat_bins, avg_fun=np.nanmean)
    return el.to(hv.Curve).clone(data=el.dframe().dropna())# * el

# colors = hv.Cycle(colorcet.b_glasbey_hv)
colors = hv.Cycle(colorcet.b_glasbey_category10)
# colors = hv.Cycle(['#282E3E', '#3F8196', '#41E2DD'])
# colors = hv.Cycle(['#2E6353', '#D1C461', '#E78A9C'])
options = [
    opts.Scatter(color=colors, tools=['hover']),
    opts.Curve(
        color=colors, line_width=5, height=350, width=500, tools=['hover'],
        xticks=[(logx, format_exponent_unicode(logx)) for logx in range(-6, -2)],
        xlim=(-6,-3.5),
        ylabel='',
    ),
    opts.RGB(),
    opts.Contours(
        show_legend=False,
        cmap=hv.Cycle([[c] for c in colors.values]),
        alpha=dim('Density').norm(), 
        line_dash='', line_width=3*np.sqrt(dim('Density').norm()),
    ),
    opts.Overlay(legend_position='bottom_left'),
]

avg = hv.NdOverlay({x: binavg(sc[x]) for x in btm_cats})

kde_args = dict(
#     levels=np.arange(0.01, 0.5, 0.02),
    bandwidth=0.3
)
kde = hv.NdOverlay({var: bivariate_kde(sc[var], **kde_args) for var in btm_cats})

l = kde * avg
l = l.opts(*options).redim(ntr0='no3_sfc_winter_uM', logstrat='strat').redim(**dims)

panel_B = l.clone().relabel('B')
panel_B
ERROR:root:Internal Python error in the inspect module.
Below is the traceback from this internal error.

Traceback (most recent call last):
  File "/Users/doppler/anaconda3/envs/nitrateflux/lib/python3.7/site-packages/IPython/core/interactiveshell.py", line 3326, in run_code
    exec(code_obj, self.user_global_ns, self.user_ns)
  File "<ipython-input-8-424344716e28>", line 19, in <module>
    btm_cats = sc.data.keys()
NameError: name 'sc' is not defined

During handling of the above exception, another exception occurred:

Traceback (most recent call last):
  File "/Users/doppler/anaconda3/envs/nitrateflux/lib/python3.7/site-packages/IPython/core/interactiveshell.py", line 2040, in showtraceback
    stb = value._render_traceback_()
AttributeError: 'NameError' object has no attribute '_render_traceback_'

During handling of the above exception, another exception occurred:

Traceback (most recent call last):
  File "/Users/doppler/anaconda3/envs/nitrateflux/lib/python3.7/site-packages/IPython/core/ultratb.py", line 1101, in get_records
    return _fixed_getinnerframes(etb, number_of_lines_of_context, tb_offset)
  File "/Users/doppler/anaconda3/envs/nitrateflux/lib/python3.7/site-packages/IPython/core/ultratb.py", line 319, in wrapped
    return f(*args, **kwargs)
  File "/Users/doppler/anaconda3/envs/nitrateflux/lib/python3.7/site-packages/IPython/core/ultratb.py", line 353, in _fixed_getinnerframes
    records = fix_frame_records_filenames(inspect.getinnerframes(etb, context))
  File "/Users/doppler/anaconda3/envs/nitrateflux/lib/python3.7/inspect.py", line 1502, in getinnerframes
    frameinfo = (tb.tb_frame,) + getframeinfo(tb, context)
  File "/Users/doppler/anaconda3/envs/nitrateflux/lib/python3.7/inspect.py", line 1460, in getframeinfo
    filename = getsourcefile(frame) or getfile(frame)
  File "/Users/doppler/anaconda3/envs/nitrateflux/lib/python3.7/inspect.py", line 696, in getsourcefile
    if getattr(getmodule(object, filename), '__loader__', None) is not None:
  File "/Users/doppler/anaconda3/envs/nitrateflux/lib/python3.7/inspect.py", line 742, in getmodule
    os.path.realpath(f)] = module.__name__
AttributeError: module has no attribute '__name__'
---------------------------------------------------------------------------

The overall statistical significance is hard to assess, but a category-wise box/whisker plot also consistently shows the ordering Basin < Shelf < Slope:

In [ ]:
hv.BoxWhisker(df, ['strat_cat', 'btm_cat'], 'ntr0').opts(xrotation=45, width=1000)

Merge panels

In [ ]:
l = panel_A + panel_B

fname = '../nb_fig/FIGURE_FN-COMP_chart_vs_NO3'

hv.save(l.opts(toolbar=None, clone=True), fname+'.html')
save_bokeh_svg_multipanel(l, fname+'.svg')
os.system(f'convert -density 200 {fname}.svg {fname}.png')